clear
%% Problem parameters
%   'EqualMinima','Himmelblau','Sixhump','ModifiedRastrigin','Vincent',
%   'Shubert','Composition','IncreasingMinima','Rastrigin','Schaffer'.
ProblemSet.FuncType = 'ModifiedRastrigin';
ProblemSet.dimension = 5;
ProblemSet.k = [2,2,2,2,2];
% ProblemSet.BaseFunc = 3;

[d, x_bound, ObjFunc, OptSet, precision_init, precision, ~, MaximalBudget, FileName]...
    = problem_setting(ProblemSet);

MaximalBudget = 1E6;
rep = 20;
iepsilon = 1;
epsilon = [0.1,0.01,0.001,0.0001];
%% problem
NewRegionNum = 2;
MaxPartitionDepth_d = ceil(log((x_bound(2,:)-x_bound(1,:))./precision(iepsilon))...
    ./log(NewRegionNum));
MaxPartitionDepth = sum(MaxPartitionDepth_d)
ClearRadius0 = (x_bound(2,:)-x_bound(1,:))./(NewRegionNum.^(MaxPartitionDepth_d));
ClearRadius = 2*min(ClearRadius0); % the maximum distance within a non-partitionable region

Problem.Dimension = d;
Problem.Domain = reshape(x_bound,[1,d*2]);
Problem.Partition = @(Region,SampleSet)Partition_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
Problem.Sampling = @(n, region)Sampling_BoxConstraint(n, region, ObjFunc);

%% optimization
% Algorithm parameters
AlgorithmM.n0 = 6;
AlgorithmM.QuantileLevel = 0.3;
AlgorithmM.NewBudget = 3;
AlgorithmM.MaxSampleSize = 15;
AlgorithmM.radius = ClearRadius;
AlgorithmM.StopCriteria = [1,MaximalBudget];
% AlgorithmM.StopCriteria = [2,size(OptSet.sol,1)];
% AlgorithmM.StopCriteria = [3,5000];

%% STAGE ONE PRS-MMO
record_gap = 1000;
times_N = zeros(floor(MaximalBudget/record_gap),rep);
for i = 1:rep
    startT = tic;
    disp(['rep: ',num2str(i)])
    [optima, RegionSet, SampleSet, RegionSampleId, time_N] ...
        = prsmmo( Problem, AlgorithmM, record_gap );
    tot_time = toc(startT);
    times_N(:,i) = time_N;
    disp(['total time used: ', num2str(tot_time)])
end

%% save data
% clear optima RegionSet SampleSet RegionSampleId time_N
% save(FileName)

%% FIGURES
times_mean = mean(times_N,2);
times_std = sqrt(var(times_N,[],2)/rep);
times_ci = times_std.*norminv(0.975);
figure()
hold on
plot(record_gap:record_gap:MaximalBudget,times_mean,'-')
% plot(record_gap:record_gap:MaximalBudget,times_mean+times_ci,'--')
% plot(record_gap:record_gap:MaximalBudget,times_mean-times_ci,'--')
hold off
xlabel('Total budget size')
ylabel('Computational time [seconds]')
set(gca,'Fontname','Times New Roman','FontSize',16);
